The data is included within this reproducible report and can be downloaded here.
# data frames with sum scales
data_scales <- data_rbt%>%
mutate(Treatment = factor(treat, levels = c("GB", "CC", "CB")),
Experitse = rowMeans(data.frame(exp_1, exp_2, exp_3,
exp_4, exp_5, exp_6), na.rm = T),
Integrity = rowMeans(data.frame(int_1, int_2, int_3, int_4), na.rm = T),
Benevolence = rowMeans(data.frame(ben_1, ben_2, ben_3, ben_4), na.rm = T),
TSM = rowMeans(data.frame(tsm_1, tsm_2, tsm_3, tsm_4), na.rm = T))
data_scales_lables <- data_rbt%>%
mutate(Treatment = factor(treat, levels = c("GB", "CC", "CB")),
Treatment = fct_recode(Treatment,
`Colored badges` = "CB",
`Control Condition` = "CC",
`Greyed out badges` = "GB"),
Experitse = rowMeans(data.frame(exp_1, exp_2, exp_3,
exp_4, exp_5, exp_6), na.rm = T),
Integrity = rowMeans(data.frame(int_1, int_2, int_3, int_4), na.rm = T),
Benevolence = rowMeans(data.frame(ben_1, ben_2, ben_3, ben_4), na.rm = T),
`Topic specific multiplism` = rowMeans(data.frame(tsm_1, tsm_2, tsm_3, tsm_4),
na.rm = T))
data_scales_wide <- data_scales%>%
select(Experitse, Integrity, Benevolence,
TSM, Treatment, session)%>%
gather(Variable, Value, -session, -Treatment)%>%
mutate(Variable2 = paste(Variable, Treatment, sep = "_"))%>%
select(-Variable, -Treatment)%>%
spread(Variable2, Value)
data_scales_lables_wide <- data_scales_lables%>%
select(Experitse, Integrity, Benevolence,
`Topic specific multiplism`, Treatment, session)%>%
gather(Variable, Value, -session, -Treatment)%>%
mutate(Variable2 = paste(Variable, Treatment, sep = "_"))%>%
select(-Variable, -Treatment)%>%
spread(Variable2, Value)
# Data of treatment check
data_tch <- data_rbt %>%
select(starts_with("tch"), meas_rep, treat)
data_tch_n <- data_tch%>%
mutate_all(function(x) ifelse(x == -999, NA, x)) %>%
filter(rowSums(is.na(data.frame(tch_1, tch_2, tch_3, tch_4, tch_5))) != 5)
PID_list <- data_rbt$session %>% unique()
length(PID_list)
## [1] 270
sum(is.na(data_rbt%>%
filter(meas_rep == 2)%>%
pull(exp_1)))
## [1] 0
data_rbt %>%
select(age, semester, sex) %>%
skim(.)
| Name | Piped data |
| Number of rows | 527 |
| Number of columns | 3 |
| _______________________ | |
| Column type frequency: | |
| numeric | 3 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| age | 21 | 0.96 | 22.89 | 2.95 | 18 | 20 | 22 | 25 | 33 | ▇▆▅▁▁ |
| semester | 21 | 0.96 | 5.86 | 3.68 | 1 | 3 | 5 | 9 | 19 | ▇▅▅▁▁ |
| sex | 21 | 0.96 | 1.71 | 0.47 | 1 | 1 | 2 | 2 | 3 | ▃▁▇▁▁ |
data_rbt$sex %>%
table(.)
## .
## 1 2 3
## 150 352 4
skewn <- function(x) DescTools::Skew(x, na.rm = T)
kurto <- function(x) DescTools::Kurt(x, na.rm = T)
maxabszvalue <- function(x) max(abs(scale(na.omit(x))))
my_skim <- skim_with(numeric = sfl(skewn, kurto, maxabszvalue))
data_scales%>%
my_skim(.)
| Name | Piped data |
| Number of rows | 527 |
| Number of columns | 34 |
| _______________________ | |
| Column type frequency: | |
| character | 2 |
| factor | 1 |
| numeric | 31 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| session | 0 | 1 | 64 | 64 | 0 | 270 | 0 |
| treat | 0 | 1 | 2 | 2 | 0 | 3 | 0 |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| Treatment | 0 | 1 | FALSE | 3 | GB: 184, CC: 173, CB: 170 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist | skewn | kurto | maxabszvalue |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| exp_1 | 0 | 1.00 | 5.43 | 1.29 | 1 | 5.00 | 6.00 | 6.0 | 7 | ▁▁▂▃▇ | -0.87 | 0.40 | 3.42 |
| exp_2 | 0 | 1.00 | 5.44 | 1.15 | 2 | 5.00 | 6.00 | 6.0 | 7 | ▂▃▆▇▅ | -0.53 | -0.21 | 2.98 |
| exp_3 | 0 | 1.00 | 5.38 | 1.26 | 2 | 5.00 | 6.00 | 6.0 | 7 | ▂▃▆▇▅ | -0.63 | -0.15 | 2.69 |
| exp_4 | 0 | 1.00 | 5.31 | 1.32 | 1 | 4.00 | 6.00 | 6.0 | 7 | ▁▁▂▃▇ | -0.74 | 0.06 | 3.25 |
| exp_5 | 0 | 1.00 | 5.10 | 1.39 | 1 | 4.00 | 5.00 | 6.0 | 7 | ▁▂▃▅▇ | -0.60 | -0.12 | 2.96 |
| exp_6 | 0 | 1.00 | 5.43 | 1.24 | 1 | 5.00 | 6.00 | 6.0 | 7 | ▁▁▂▃▇ | -0.68 | 0.02 | 3.57 |
| int_1 | 0 | 1.00 | 5.32 | 1.24 | 1 | 4.00 | 6.00 | 6.0 | 7 | ▁▁▃▃▇ | -0.49 | -0.35 | 3.48 |
| int_2 | 0 | 1.00 | 5.33 | 1.31 | 1 | 4.00 | 6.00 | 6.0 | 7 | ▁▁▃▃▇ | -0.47 | -0.50 | 3.30 |
| int_3 | 0 | 1.00 | 5.23 | 1.17 | 1 | 4.00 | 5.00 | 6.0 | 7 | ▁▁▅▅▇ | -0.10 | -0.78 | 3.60 |
| int_4 | 0 | 1.00 | 5.32 | 1.22 | 1 | 4.00 | 6.00 | 6.0 | 7 | ▁▁▃▃▇ | -0.45 | -0.37 | 3.53 |
| ben_1 | 0 | 1.00 | 5.23 | 1.21 | 1 | 4.00 | 5.00 | 6.0 | 7 | ▁▁▅▅▇ | -0.23 | -0.45 | 3.51 |
| ben_2 | 0 | 1.00 | 5.17 | 1.26 | 1 | 4.00 | 5.00 | 6.0 | 7 | ▁▁▆▅▇ | -0.22 | -0.34 | 3.31 |
| ben_3 | 0 | 1.00 | 5.27 | 1.26 | 1 | 4.00 | 5.00 | 6.0 | 7 | ▁▁▃▃▇ | -0.63 | 0.26 | 3.37 |
| ben_4 | 0 | 1.00 | 5.00 | 1.23 | 1 | 4.00 | 5.00 | 6.0 | 7 | ▁▂▆▅▇ | -0.09 | -0.56 | 3.26 |
| meas_rep | 0 | 1.00 | 1.49 | 0.50 | 1 | 1.00 | 1.00 | 2.0 | 2 | ▇▁▁▁▇ | 0.05 | -2.00 | 1.02 |
| tsm_1 | 0 | 1.00 | 1.87 | 0.87 | 1 | 1.00 | 2.00 | 2.0 | 4 | ▇▇▁▃▁ | 0.69 | -0.34 | 2.46 |
| tsm_2 | 0 | 1.00 | 2.05 | 0.95 | 1 | 1.00 | 2.00 | 3.0 | 4 | ▇▇▁▅▂ | 0.51 | -0.73 | 2.06 |
| tsm_3 | 0 | 1.00 | 2.05 | 0.97 | 1 | 1.00 | 2.00 | 3.0 | 4 | ▇▇▁▃▂ | 0.62 | -0.62 | 2.01 |
| tsm_4 | 0 | 1.00 | 2.29 | 1.01 | 1 | 1.00 | 2.00 | 3.0 | 4 | ▆▇▁▆▃ | 0.23 | -1.06 | 1.69 |
| age | 21 | 0.96 | 22.89 | 2.95 | 18 | 20.00 | 22.00 | 25.0 | 33 | ▇▆▅▁▁ | 0.72 | -0.06 | 3.43 |
| semester | 21 | 0.96 | 5.86 | 3.68 | 1 | 3.00 | 5.00 | 9.0 | 19 | ▇▅▅▁▁ | 0.57 | -0.44 | 3.57 |
| sex | 21 | 0.96 | 1.71 | 0.47 | 1 | 1.00 | 2.00 | 2.0 | 3 | ▃▁▇▁▁ | -0.70 | -0.91 | 2.74 |
| tch_1 | 0 | 1.00 | -107.78 | 313.71 | -999 | 1.00 | 2.00 | 4.0 | 4 | ▁▁▁▁▇ | -2.48 | 4.18 | 2.84 |
| tch_2 | 0 | 1.00 | -265.55 | 443.71 | -999 | -999.00 | 1.00 | 3.5 | 4 | ▃▁▁▁▇ | -1.05 | -0.91 | 1.65 |
| tch_3 | 0 | 1.00 | -278.85 | 450.45 | -999 | -999.00 | 1.00 | 4.0 | 4 | ▃▁▁▁▇ | -0.97 | -1.06 | 1.60 |
| tch_4 | 0 | 1.00 | -111.51 | 318.42 | -999 | 1.00 | 2.00 | 4.0 | 4 | ▁▁▁▁▇ | -2.42 | 3.89 | 2.79 |
| tch_5 | 0 | 1.00 | -276.95 | 449.52 | -999 | -999.00 | 1.00 | 3.0 | 4 | ▃▁▁▁▇ | -0.98 | -1.04 | 1.61 |
| Experitse | 0 | 1.00 | 5.35 | 1.12 | 2 | 4.67 | 5.50 | 6.0 | 7 | ▁▂▅▇▅ | -0.59 | -0.09 | 3.00 |
| Integrity | 0 | 1.00 | 5.30 | 1.08 | 2 | 4.50 | 5.50 | 6.0 | 7 | ▁▃▆▇▅ | -0.28 | -0.60 | 3.05 |
| Benevolence | 0 | 1.00 | 5.17 | 1.06 | 1 | 4.25 | 5.25 | 6.0 | 7 | ▁▁▆▇▆ | -0.16 | -0.16 | 3.93 |
| TSM | 0 | 1.00 | 2.06 | 0.67 | 1 | 1.50 | 2.00 | 2.5 | 4 | ▆▆▇▂▁ | 0.33 | -0.25 | 2.88 |
First we specified two consecutive threedimensional CFA models
# onedimensional model
cfa_meti_model_1d <- "exp =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6 +
int_1 + int_2 + int_3 + int_4 +
ben_1 + ben_2 + ben_3 + ben_4"
cfa1d_meti_1 <- cfa(cfa_meti_model_1d, data = data_rbt%>%filter(meas_rep == 1))
cfa1d_meti_2 <- cfa(cfa_meti_model_1d, data = data_rbt%>%filter(meas_rep == 2))
cfa_meti_model_1 <- "exp =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6
int =~ int_1 + int_2 + int_3 + int_4
ben =~ ben_1 + ben_2 + ben_3 + ben_4
int_1 ~~ int_2"
cfa_meti_model_2 <- "exp =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6
int =~ int_1 + int_2 + int_3 + int_4
ben =~ ben_1 + ben_2 + ben_3 + ben_4
ben_1 ~~ ben_3
ben_1 ~~ ben_2"
cfa_meti_1 <- cfa(cfa_meti_model_1, data = data_rbt%>%filter(meas_rep == 1))
# define a function which prints the fit
fpf <- function(x){
fm_tmp <- fitmeasures(x)
return(cat(sprintf(
"χ^2^ = %s, _df_ = %s, CFI = %s, TLI = %s, RMSEA = %s, SRMR = %s, SRMR~between~ = %s, SRMR~within~ = %s",
round(fm_tmp[c("chisq")],3),
fm_tmp[c("df")],
round(fm_tmp[c("cfi")],3),
round(fm_tmp[c("tli")],3),
round(fm_tmp[c("rmsea")],3),
round(fm_tmp[c("srmr")],3),
round(fm_tmp[c("srmr_between")],3),
round(fm_tmp[c("srmr_within")],3))))
}
# print the fit for cfa_meti_1
fpf(cfa1d_meti_1)
χ2 = 588.932, df = 77, CFI = 0.811, TLI = 0.776, RMSEA = 0.157, SRMR = 0.084, SRMRbetween = NA, SRMRwithin = NA
fpf(cfa_meti_1)
χ2 = 194.174, df = 73, CFI = 0.955, TLI = 0.944, RMSEA = 0.078, SRMR = 0.049, SRMRbetween = NA, SRMRwithin = NA
cfa_meti_2 <- cfa(cfa_meti_model_2, data = data_rbt%>%filter(meas_rep == 2))
fpf(cfa1d_meti_2)
χ2 = 936.555, df = 77, CFI = 0.759, TLI = 0.715, RMSEA = 0.208, SRMR = 0.099, SRMRbetween = NA, SRMRwithin = NA
fpf(cfa_meti_2)
χ2 = 212.262, df = 72, CFI = 0.961, TLI = 0.95, RMSEA = 0.087, SRMR = 0.04, SRMRbetween = NA, SRMRwithin = NA
anova(cfa1d_meti_1, cfa_meti_1)
## Chi-Squared Difference Test
##
## Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
## cfa_meti_1 73 9738.4 9853.5 194.17
## cfa1d_meti_1 77 10125.1 10225.9 588.93 394.76 4 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(cfa1d_meti_2, cfa_meti_2)
## Chi-Squared Difference Test
##
## Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
## cfa_meti_2 72 8522.8 8639.9 212.26
## cfa1d_meti_2 77 9237.1 9336.5 936.55 724.29 5 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
In an next step we ran a two-level CFA …
# onedimensional model
mcfa1d_meti_model <- "level: 1
meti =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6 +
int_1 + int_2 + int_3 + int_4 +
ben_1 + ben_2 + ben_3 + ben_4
level: 2
meti =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6 +
int_1 + int_2 + int_3 + int_4 +
ben_1 + ben_2 + ben_3 + ben_4"
mcfa_meti_model <- "level: 1
exp_w =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6
int_w =~ int_1 + int_2 + int_3 + int_4
ben_w =~ ben_1 + ben_2 + ben_3 + ben_4
int_1 ~~ int_2
int_3 ~~ ben_4
ben_1 ~~ ben_2
level: 2
exp_b =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6
int_b =~ int_1 + int_2 + int_3 + int_4
ben_b =~ ben_1 + ben_2 + ben_3 + ben_4"
mcfa1d_meti <- cfa(mcfa1d_meti_model, data = data_rbt, cluster = "session")
mcfa_meti <- cfa(mcfa_meti_model, data = data_rbt, cluster = "session")
fpf(mcfa1d_meti)
χ2 = 764.777, df = 154, CFI = 0.894, TLI = 0.875, RMSEA = 0.087, SRMR = 0.271, SRMRbetween = 0.146, SRMRwithin = 0.125
fpf(mcfa_meti)
χ2 = 277.759, df = 145, CFI = 0.977, TLI = 0.971, RMSEA = 0.042, SRMR = 0.172, SRMRbetween = 0.091, SRMRwithin = 0.081
anova(mcfa1d_meti, mcfa_meti)
## Chi-Squared Difference Test
##
## Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
## mcfa_meti 145 18147 18484 277.76
## mcfa1d_meti 154 18616 18915 764.78 487.02 9 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# First Measurement
## Expertise
MBESS::ci.reliability(data_rbt%>%filter(meas_rep == 1) %>% select(starts_with("exp")))$est
## [1] 0.9247211
# Integrity
MBESS::ci.reliability(data_rbt%>%filter(meas_rep == 1) %>% select(starts_with("int")))$est
## [1] 0.8806539
# Benevolence
MBESS::ci.reliability(data_rbt%>%filter(meas_rep == 1) %>% select(starts_with("ben")))$est
## [1] 0.8513474
# Second Measurement
## Expertise
MBESS::ci.reliability(data_rbt%>%filter(meas_rep == 2) %>% select(starts_with("exp")))$est
## [1] 0.9490452
# Integrity
MBESS::ci.reliability(data_rbt%>%filter(meas_rep == 2) %>% select(starts_with("int")))$est
## [1] 0.9120384
# Benevolence
MBESS::ci.reliability(data_rbt%>%filter(meas_rep == 2) %>% select(starts_with("ben")))$est
## [1] 0.9003258
First we specified two consecutive onedimensional CFA models
cfa_tsm_model <- "tsm =~ a*tsm_1 + a*tsm_2 + a*tsm_3 + a*tsm_4
tsm_2 ~~ tsm_4"
cfa_tsm_1 <- cfa(cfa_tsm_model, data = data_rbt%>%filter(meas_rep == 1))
fpf(cfa_tsm_1)
χ2 = 5.302, df = 4, CFI = 0.991, TLI = 0.987, RMSEA = 0.035, SRMR = 0.048, SRMRbetween = NA, SRMRwithin = NA
cfa_tsm_2 <- cfa(cfa_tsm_model, data = data_rbt%>%filter(meas_rep == 2))
fpf(cfa_tsm_2)
χ2 = 6.072, df = 4, CFI = 0.99, TLI = 0.985, RMSEA = 0.045, SRMR = 0.055, SRMRbetween = NA, SRMRwithin = NA
In an next step, we ran a two-level CFA …
mcfa_tsm_model <- "level: 1
tsm_w =~ tsm_1 + tsm_2 + tsm_3 + tsm_4
tsm_2 ~~ tsm_4
level: 2
tsm_b =~ tsm_1 + tsm_2 + tsm_3 + tsm_4
tsm_2 ~~ 0*tsm_2"
mcfa_tsm <- cfa(mcfa_tsm_model, data = data_rbt, cluster = "session")
fpf(mcfa_tsm)
χ2 = 4.422, df = 4, CFI = 0.999, TLI = 0.996, RMSEA = 0.014, SRMR = 0.109, SRMRbetween = 0.09, SRMRwithin = 0.019
# First Measurement
MBESS::ci.reliability(data_rbt %>% filter(meas_rep ==1) %>% select(starts_with("tsm")))$est
## [1] 0.6509044
# Second Measurement
MBESS::ci.reliability(data_rbt %>% filter(meas_rep ==2) %>% select(starts_with("tsm")))$est
## [1] 0.7125262
We specified two consecutive onedimensional CFA models
cfa_tch_model <- "tch =~ tch_1 + tch_2 + tch_3 + tch_4 + tch_5
tch_1 ~~ tch_4
tch_3 ~~ tch_5"
cfa_tch_1 <- cfa(cfa_tch_model,
data = data_tch_n%>%
filter(meas_rep == 1))
fpf(cfa_tch_1)
χ2 = 2.412, df = 3, CFI = 1, TLI = 1.003, RMSEA = 0, SRMR = 0.006, SRMRbetween = NA, SRMRwithin = NA
cfa_tch_2 <- cfa(cfa_tch_model, data = data_tch_n%>%filter(meas_rep == 2))
fpf(cfa_tch_2)
χ2 = 6.538, df = 3, CFI = 0.997, TLI = 0.991, RMSEA = 0.083, SRMR = 0.005, SRMRbetween = NA, SRMRwithin = NA
# First Measurement
MBESS::ci.reliability(data_rbt %>% filter(meas_rep == 2) %>% select(starts_with("tch")))$est
## [1] 0.8625644
# Second Measurement
MBESS::ci.reliability(data_rbt %>% filter(meas_rep == 2) %>% select(starts_with("tch")))$est
## [1] 0.8625644
tibble(`1d CFA METI 1` = fitmeasures(cfa1d_meti_1)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`1d CFA METI 2` = fitmeasures(cfa1d_meti_2)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`3d CFA METI 1` = fitmeasures(cfa_meti_1)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`3d CFA METI 2` = fitmeasures(cfa_meti_2)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`1d MCFA METI` = fitmeasures(mcfa1d_meti)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`3d MCFA METI` = fitmeasures(mcfa_meti)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`1d CFA TSM 1` = fitmeasures(cfa_tsm_1)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`1d CFA TSM 2` = fitmeasures(cfa_tsm_2)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`1d MCFA TSM` = fitmeasures(mcfa_tsm)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`1d CFA TCH 1` = fitmeasures(cfa_tch_1)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`1d CFA TCH 2` = fitmeasures(cfa_tch_2)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
rownames = c("χ^2^", "_df_", "CFI", "TLI", "RMSEA", "SRMR", "SRMR~between~", "SRMR~within~", "BIC", "AIC"))%>%
column_to_rownames(var = "rownames")%>%
knitr::kable(., digits = 3)
| 1d CFA METI 1 | 1d CFA METI 2 | 3d CFA METI 1 | 3d CFA METI 2 | 1d MCFA METI | 3d MCFA METI | 1d CFA TSM 1 | 1d CFA TSM 2 | 1d MCFA TSM | 1d CFA TCH 1 | 1d CFA TCH 2 | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| χ2 | 588.932 | 936.555 | 194.174 | 212.262 | 764.777 | 277.759 | 5.302 | 6.072 | 4.422 | 2.412 | 6.538 |
| df | 77.000 | 77.000 | 73.000 | 72.000 | 154.000 | 145.000 | 4.000 | 4.000 | 4.000 | 3.000 | 3.000 |
| CFI | 0.811 | 0.759 | 0.955 | 0.961 | 0.894 | 0.977 | 0.991 | 0.990 | 0.999 | 1.000 | 0.997 |
| TLI | 0.776 | 0.715 | 0.944 | 0.950 | 0.875 | 0.971 | 0.987 | 0.985 | 0.996 | 1.003 | 0.991 |
| RMSEA | 0.157 | 0.208 | 0.078 | 0.087 | 0.087 | 0.042 | 0.035 | 0.045 | 0.014 | 0.000 | 0.083 |
| SRMR | 0.084 | 0.099 | 0.049 | 0.040 | 0.271 | 0.172 | 0.048 | 0.055 | 0.109 | 0.006 | 0.005 |
| SRMRbetween | NA | NA | NA | NA | 0.146 | 0.091 | NA | NA | 0.090 | NA | NA |
| SRMRwithin | NA | NA | NA | NA | 0.125 | 0.081 | NA | NA | 0.019 | NA | NA |
| BIC | 10225.876 | 9336.494 | 9853.512 | 8639.946 | 18914.759 | 18484.147 | 2791.127 | 2653.377 | 5445.587 | 1547.731 | 1750.473 |
| AIC | 10125.120 | 9237.119 | 9738.362 | 8522.826 | 18616.055 | 18147.038 | 2769.537 | 2632.082 | 5360.243 | 1512.868 | 1712.703 |
res_tch_data <- data_tch%>%
gather(variable, value, starts_with("tch_"))%>%
group_by(treat, variable, value)%>%
mutate(value = as.character(value)) %>%
summarize(freq = n())%>%
ungroup()%>%
mutate(treat = case_when(treat == "GB" ~ "Greyed out badges",
treat == "CB" ~ "Colored badges",
treat == "CC" ~ "Control condition",
T ~ treat),
value = case_when(value =="-999" ~ "don't know",
T ~ value),
variable = case_when(variable == "tch_1" ~ "item 1",
variable == "tch_2" ~ "item 2",
variable == "tch_3" ~ "item 3",
variable == "tch_4" ~ "item 4",
variable == "tch_5" ~ "item 5"),
Frequency = freq)
## `summarise()` regrouping output by 'treat', 'variable' (override with `.groups` argument)
res_tch_plot <- ggplot(res_tch_data, aes(variable, value)) +
geom_point(aes(size = Frequency, color = Frequency), shape = 15) +
scale_size_continuous(range = c(3,15)) +
scale_color_gradient(low = "grey95", high = "grey65") +
guides(color=guide_legend(), size = guide_legend()) +
facet_wrap(~treat) +
theme_ipsum_ps() +
ggtitle("Results of the treatment check", "Frequency per item and experimental condition") +
ylab("") +
xlab("")
res_tch_plot
res_tch_plot_pub <- ggplot(res_tch_data %>%
mutate(treat = case_when(treat == "Greyed out badges" ~ "Grey badges",
TRUE ~ treat)),
aes(variable, value)) +
geom_point(aes(size = Frequency, color = Frequency), shape = 15) +
scale_size_continuous(range = c(3,15)) +
scale_color_gradient(low = "grey95", high = "grey65") +
guides(color=guide_legend(), size = guide_legend()) +
facet_wrap(~treat) +
theme_ipsum_ps() +
ylab("") +
xlab("")
#ggsave("8_reproducible_documentation_of_analyses/Fig/Fig2.jpg", res_tch_plot_pub, width = 130, height = 50, units = "mm", dpi = 300, scale = 2.4)
res_tch_data_A <- data_tch_n%>%
filter(treat != "CC")%>%
filter(tch_1 != -999)
effsize::VD.A(tch_1 ~ treat, data = res_tch_data_A)
##
## Vargha and Delaney A
##
## A estimate: 0.837694 (large)
data_scales_lables%>%
gather(Variable, Value, Experitse, Integrity, Benevolence)%>%
ggplot(., aes(x = Treatment, y = Value)) +
geom_violin(adjust = 1.5) +
stat_summary(fun.data = "mean_sdl", fun.args = list(mult = 1)) +
coord_flip() +
facet_wrap(~ Variable, nrow = 1) +
labs(title = "Graphical overview (Hyp 1)",
subtitle = "Violinplots and means ± 1*SD",
caption = "") +
ylim(1,7) +
hrbrthemes::theme_ipsum_ps()
# Export data for combined plot
#data_scales_lables %>%
# mutate(Study = "Study 1") %>%
# write_csv(., "8_reproducible_documentation_of_analyses/Fig/data_scales_lables_study_1.csv")
A_GB_CC <- data_scales_lables%>%
filter(treat != "CB")%>%
mutate(treat = as.character(treat))%>%
effsize::VD.A(Integrity ~ treat, data = .)%>%
.$estimate
A_CC_CB <- data_scales_lables%>%
filter(treat != "GB")%>%
mutate(treat = as.character(treat))%>%
effsize::VD.A(Integrity ~ treat, data = .)%>%
.$estimate
A_GB_CB <- data_scales_lables%>%
filter(treat != "CC")%>%
mutate(treat = as.character(treat))%>%
effsize::VD.A(Integrity ~ treat, data = .)%>%
.$estimate
| GB | CC | |
|---|---|---|
| CC | A = 0.59, d = 0.32 | |
| CB | A = 0.66, d = 0.57 | A = 0.58, d = 0.29 |
plot_hyp2_1 <- ggplot(data_scales_lables,
aes(x=`Topic specific multiplism`, y = Integrity)) +
geom_jitter() +
facet_wrap(~ Treatment, nrow = 1) +
labs(title = "Graphical overview (Hyp 2/3)",
subtitle = "Jitter plot per treatment") +
hrbrthemes::theme_ipsum_ps()
plot_hyp2_1 + stat_smooth()
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
plot_hyp2_1 + stat_smooth(method = "lm")
## `geom_smooth()` using formula 'y ~ x'
Spearman and Kendall correlations:
r_GB <- round(cor(data_scales_lables%>%
filter(treat == "GB")%>%
pull(Integrity),
data_scales_lables%>%
filter(treat == "GB")%>%
pull(`Topic specific multiplism`),
method = "pearson", use = "pairwise.complete.obs"), 2)
t_GB <- round(cor(data_scales_lables%>%
filter(treat == "GB")%>%
pull(Integrity),
data_scales_lables%>%
filter(treat == "GB")%>%
pull(`Topic specific multiplism`),
method = "kendall", use = "pairwise.complete.obs"), 2)
r_CC <- round(cor(data_scales_lables%>%
filter(treat == "CC")%>%
pull(Integrity),
data_scales_lables%>%
filter(treat == "CC")%>%
pull(`Topic specific multiplism`),
method = "pearson", use = "pairwise.complete.obs"), 2)
t_CC <- round(cor(data_scales_lables%>%
filter(treat == "CC")%>%
pull(Integrity),
data_scales_lables%>%
filter(treat == "CC")%>%
pull(`Topic specific multiplism`),
method = "kendall", use = "pairwise.complete.obs"), 2)
r_CB <- round(cor(data_scales_lables%>%
filter(treat == "CB")%>%
pull(Integrity),
data_scales_lables%>%
filter(treat == "CB")%>%
pull(`Topic specific multiplism`),
method = "pearson", use = "pairwise.complete.obs"), 2)
t_CB <- round(cor(data_scales_lables%>%
filter(treat == "CB")%>%
pull(Integrity),
data_scales_lables%>%
filter(treat == "CB")%>%
pull(`Topic specific multiplism`),
method = "kendall", use = "pairwise.complete.obs"), 2)
| GB | CC | CB | |
|---|---|---|---|
| \(r(integrity, multiplism)\) | -0.46 | -0.36 | -0.35 |
| \(\tau(integrity, multiplism)\) | -0.33 | -0.24 | -0.26 |
data_scales_lables%>%
mutate(Treatment = case_when(treat == "GB" ~ "Greyed out badges",
treat == "CB" ~ "Colored badges",
treat == "CC" ~ "Control condition",
T ~ "treat"))%>%
ggplot(., aes(x = Treatment, y = `Topic specific multiplism`)) +
geom_violin(adjust = 1.5, alpha = .5) +
stat_summary(fun.data = "mean_sdl", fun.args = list(mult = 1)) +
coord_flip() +
labs(title = "Graphical overview (Hyp 4)",
subtitle = "Violinplots and means ± 1*SD",
caption = "") +
xlab("") +
ylim(1,4) +
hrbrthemes::theme_ipsum_ps()
fig4 <- data_scales_lables%>%
mutate(Treatment = case_when(treat == "GB" ~ "Greyed out badges",
treat == "CB" ~ "Colored badges",
treat == "CC" ~ "Control condition",
T ~ "treat"))%>%
ggplot(., aes(x = Treatment, y = `Topic specific multiplism`)) +
geom_violin(adjust = 1.5, alpha = .5) +
stat_summary(fun.data = "mean_sdl", fun.args = list(mult = 1)) +
coord_flip() +
xlab("") +
ylim(1,4) +
hrbrthemes::theme_ipsum_ps()
A_mult_GB_CC <- effsize::VD.A(`Topic specific multiplism` ~ treat,
data = data_scales_lables%>%
filter(treat != "CB")%>%
mutate(treat = as.character(treat)))$estimate
d_mult_GB_CC <- qnorm(A_mult_GB_CC)*sqrt(2)
A_mult_CC_CB <- effsize::VD.A(`Topic specific multiplism` ~ treat,
data = data_scales_lables%>%
filter(treat != "GB")%>%
mutate(treat = as.character(treat)))$estimate
d_mult_CC_CB <- qnorm(A_mult_CC_CB)*sqrt(2)
A_mult_GB_CB <- effsize::VD.A(`Topic specific multiplism` ~ treat,
data = data_scales_lables%>%
filter(treat != "CC")%>%
mutate(treat = as.character(treat)))$estimate
d_mult_GB_CB <- qnorm(A_mult_GB_CB)*sqrt(2)
| GB | CC | |
|---|---|---|
| CC | A = 0.43, d = -0.26 | |
| CB | A = 0.43, d = -0.25 | A = 0.5, d = 0.01 |
All of the following analyses are based on the data frame object data_scales_wide which is why we describe it here somewhat more detailed.
All analyses are based on measured variables Integrity (Integrity, is information source sincere, honest, just, unselfish and fair?) and Topic Specific Multiplism (TSM, is knowledge and knowing about a topic arbitrary?). As this data set is in wide format, the experimental conditions are encoded wtihin the variable names:
GB means, that the participants of the study could learn from the grey out badges (ans corresponding explanations) within the abstract, that the authors of the study denied to use Open PracticesCC means, that the participants of the study could not learn if or if not the authors of the study used Open PracticesCBmeans, that the participants of the study could learn from the colored badges (ans corresponding explanations) within the abstract, that the authors of the study used Open PracticesFinally, the variable session identified the study participants.
If we look descriptively at these variables:
data_scales_wide%>%
my_skim(.)
| Name | Piped data |
| Number of rows | 270 |
| Number of columns | 13 |
| _______________________ | |
| Column type frequency: | |
| character | 1 |
| numeric | 12 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| session | 0 | 1 | 64 | 64 | 0 | 270 | 0 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist | skewn | kurto | maxabszvalue |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Benevolence_CB | 100 | 0.63 | 5.44 | 0.99 | 2.75 | 4.75 | 5.25 | 6.25 | 7 | ▁▂▇▅▆ | -0.08 | -0.71 | 2.71 |
| Benevolence_CC | 97 | 0.64 | 5.20 | 0.96 | 2.75 | 4.50 | 5.25 | 6.00 | 7 | ▁▅▇▆▃ | 0.03 | -0.61 | 2.55 |
| Benevolence_GB | 86 | 0.68 | 4.89 | 1.15 | 1.00 | 4.00 | 4.75 | 5.75 | 7 | ▁▁▇▆▅ | -0.14 | -0.02 | 3.39 |
| Experitse_CB | 100 | 0.63 | 5.63 | 1.10 | 2.33 | 4.88 | 5.83 | 6.50 | 7 | ▁▂▃▇▇ | -0.78 | 0.13 | 2.99 |
| Experitse_CC | 97 | 0.64 | 5.35 | 0.97 | 2.50 | 4.67 | 5.33 | 6.00 | 7 | ▁▃▆▇▅ | -0.25 | -0.48 | 2.94 |
| Experitse_GB | 86 | 0.68 | 5.09 | 1.19 | 2.00 | 4.29 | 5.33 | 6.00 | 7 | ▂▂▅▇▃ | -0.59 | -0.32 | 2.59 |
| Integrity_CB | 100 | 0.63 | 5.62 | 0.98 | 3.25 | 5.00 | 5.75 | 6.50 | 7 | ▂▅▇▇▇ | -0.31 | -0.83 | 2.42 |
| Integrity_CC | 97 | 0.64 | 5.34 | 0.97 | 2.75 | 4.50 | 5.50 | 6.00 | 7 | ▁▃▇▇▅ | -0.22 | -0.62 | 2.66 |
| Integrity_GB | 86 | 0.68 | 4.97 | 1.18 | 2.00 | 4.00 | 5.00 | 6.00 | 7 | ▂▅▆▇▃ | -0.09 | -0.76 | 2.53 |
| TSM_CB | 100 | 0.63 | 2.01 | 0.67 | 1.00 | 1.50 | 2.00 | 2.50 | 4 | ▇▇▇▂▁ | 0.38 | -0.34 | 2.95 |
| TSM_CC | 97 | 0.64 | 1.99 | 0.61 | 1.00 | 1.50 | 2.00 | 2.50 | 4 | ▆▆▇▁▁ | 0.25 | -0.12 | 3.30 |
| TSM_GB | 86 | 0.68 | 2.18 | 0.72 | 1.00 | 1.75 | 2.25 | 2.56 | 4 | ▅▆▇▂▂ | 0.24 | -0.46 | 2.54 |
library(naniar)
##
## Attaching package: 'naniar'
## The following object is masked from 'package:skimr':
##
## n_complete
visdat::vis_miss(data_scales_wide%>%select(-session)) +
ggtitle("Missingness per Variable") +
theme(plot.margin = margin(1, 2, 1, 1, "cm"))
Integrity_GBlibrary(patchwork)
ggplot(data_scales_wide, aes(Integrity_GB, Integrity_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, Integrity_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, Integrity_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB, Integrity_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC, Integrity_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB, Integrity_GB)) +
geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
theme_ipsum_ps()
Integrity_CCggplot(data_scales_wide, aes(Integrity_GB, Integrity_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, Integrity_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, Integrity_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB, Integrity_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC, Integrity_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB, Integrity_CC)) +
geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
theme_ipsum_ps()
Integrity_CBggplot(data_scales_wide, aes(Integrity_GB, Integrity_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, Integrity_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, Integrity_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB, Integrity_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC, Integrity_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB, Integrity_CB)) +
geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
theme_ipsum_ps()
TSM_GBggplot(data_scales_wide, aes(Integrity_GB, TSM_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, TSM_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, TSM_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB, TSM_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC, TSM_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB, TSM_GB)) +
geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
theme_ipsum_ps()
TSM_CCggplot(data_scales_wide, aes(Integrity_GB, TSM_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, TSM_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, TSM_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB, TSM_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC, TSM_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB, TSM_CC)) +
geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
theme_ipsum_ps()
TSM_CBggplot(data_scales_wide, aes(Integrity_GB, TSM_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, TSM_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, TSM_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB, TSM_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC, TSM_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB, TSM_CB)) +
geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
theme_ipsum_ps()
d_matrix_missings <- miceadds::mi_dstat(data_scales_wide%>%select(starts_with("Int"), starts_with("TSM")))%>%
round(.,4)
knitr::kable(d_matrix_missings, caption = "Boxplot of Cohen's d of missing/not-missing per variable")
| Integrity_CB | Integrity_CC | Integrity_GB | TSM_CB | TSM_CC | TSM_GB | |
|---|---|---|---|---|---|---|
| Integrity_CB | NaN | -0.2103 | -0.0096 | NaN | -0.1009 | -0.3312 |
| Integrity_CC | -0.1596 | NaN | -0.0074 | 0.1056 | NaN | 0.3622 |
| Integrity_GB | 0.1045 | 0.1696 | NaN | -0.1339 | 0.1018 | NaN |
| TSM_CB | NaN | -0.2103 | -0.0096 | NaN | -0.1009 | -0.3312 |
| TSM_CC | -0.1596 | NaN | -0.0074 | 0.1056 | NaN | 0.3622 |
| TSM_GB | 0.1045 | 0.1696 | NaN | -0.1339 | 0.1018 | NaN |
boxplot(as.vector(d_matrix_missings))
title("Boxplot of Cohen's d of missing/not-missing per variable")
M <- 1000
out <- mice(data = data_scales_wide%>%select(-session),
m = M,
meth=c("norm","norm","norm",
"norm","norm","norm",
"norm","norm","norm",
"norm","norm","norm"),
diagnostics = FALSE,
printFlag = F,
seed = 83851)
out_first10 <- mice(data = data_scales_wide%>%select(-session),
m = 10,
meth=c("norm","norm","norm",
"norm","norm","norm",
"norm","norm","norm",
"norm","norm","norm"),
diagnostics = FALSE,
printFlag = F,
seed = 83851)
densityplot(out_first10)
plot(out_first10)
# Set up the matrices for the estimates ##############
# setup of matrices to store multiple estimates
mulest_hyp1 <- matrix(0,nrow=M,ncol=3)
# and covariance matrices
covwithin_hyp1 <- matrix(0,nrow=3,ncol=3)
# Estimate the coefficients for each data frame ######
for(i in 1:M) {
within_hyp1 <- lm(cbind(Integrity_GB,Integrity_CC,Integrity_CB)~1,
data=mice::complete(out,i)) # estimate the means of the three variables
mulest_hyp1[i,]<-coef(within_hyp1)[1:3] # store these means in the matrix `mulres`
covwithin_hyp1<-covwithin_hyp1 + 1/M * vcov(within_hyp1)[1:3,1:3] # compute the averages
}
# Compute the average of the estimates ###############
estimates_hyp1 <- colMeans(mulest_hyp1)
names(estimates_hyp1) <- c("Integrity_GB","Integrity_CC","Integrity_CB")
covbetween_hyp1 <- cov(mulest_hyp1) # between covariance matrix
covariance_hyp1 <- covwithin_hyp1 + (1+1/M)*covbetween_hyp1 # total variance
# Determine the effective and real sample sizes ######
samp_hyp1 <- nrow(data_scales_wide) # real sample size
nucom_hyp1<-samp_hyp1-length(estimates_hyp1)
# corresponds to Equation (X) in Hoijtink, Gu, Mulder, & Rosseel (2019) ...
lam_hyp1 <- (1+1/M)*(1/length(estimates_hyp1))*
sum(diag(covbetween_hyp1 %*%
MASS::ginv(covariance_hyp1))) # ... (43)
nuold_hyp1<-(M-1)/(lam_hyp1^2) # ... (44)
nuobs_hyp1<-(nucom_hyp1+1)/(nucom_hyp1+3)*nucom_hyp1*(1-lam_hyp1) # ... (46)
nu_hyp1<- nuold_hyp1*nuobs_hyp1/(nuold_hyp1+nuobs_hyp1) # ... (47)
fracmis_hyp1 <- (nu_hyp1+1)/(nu_hyp1+3)*lam_hyp1 + 2/(nu_hyp1+3) # ... (48)
neff_hyp1<-samp_hyp1-samp_hyp1*fracmis_hyp1 # = 172 approx. 2/3* 270
# coerce `covariance` to a list
covariance_hyp1<-list(covariance_hyp1)
# Test the hypotheses with bain ######################
set.seed(6346)
results_hyp1 <- bain(estimates_hyp1,
"Integrity_GB<Integrity_CC<Integrity_CB;
Integrity_GB=Integrity_CC=Integrity_CB;
Integrity_GB<Integrity_CC=Integrity_CB",
n = neff_hyp1, Sigma=covariance_hyp1,
group_parameters=3, joint_parameters = 0)
print(results_hyp1)
## Bayesian informative hypothesis testing for an object of class numeric:
##
## Fit Com BF.u BF.c PMPa PMPb
## H1 0.999 0.177 5.656 3878.238 0.970 0.828
## H2 0.000 0.279 0.000 0.000 0.000 0.000
## H3 0.045 0.257 0.174 0.174 0.030 0.025
## Hu 0.146
##
## Hypotheses:
## H1: Integrity_GB<Integrity_CC<Integrity_CB
## H2: Integrity_GB=Integrity_CC=Integrity_CB
## H3: Integrity_GB<Integrity_CC=Integrity_CB
##
## Note: BF.u denotes the Bayes factor of the hypothesis at hand versus the unconstrained hypothesis Hu. BF.c denotes the Bayes factor of the hypothesis at hand versus its complement.
summary(results_hyp1)
## Parameter n Estimate lb ub
## 1 Integrity_GB 168.9123 4.992512 4.830189 5.154835
## 2 Integrity_CC 168.9123 5.327959 5.189725 5.466193
## 3 Integrity_CB 168.9123 5.585524 5.443602 5.727446
results_hyp1$BFmatrix
## H1 H2 H3
## H1 1.000000e+00 61491369 3.254405e+01
## H2 1.626244e-08 1 5.292458e-07
## H3 3.072758e-02 1889481 1.000000e+00
path_mod <- "Integrity_GB ~ TSM_GB
Integrity_CC ~ TSM_CC
Integrity_CB ~ TSM_CB"
# Set up the matrices for the estimates ##############
best_hyp2 <- matrix(0,nrow=M,ncol=3) # setup of matrices to store multiple estimates
covwithin_hyp2 <- matrix(0,nrow=3,ncol=3) # and covariance matrices
# Estimate the coefficients for each data frame ######
for(i in 1:M) {
path_fit <- sem(path_mod,
data = mice::complete(out, i),
std.lv = T) # estimate the path coefficients
best_hyp2[i,] <- parameterestimates(path_fit, standardized = T)%>% # store path coefficients
filter(op == "~")%>%
pull(std.all)
covwithin_hyp2 <- covwithin_hyp2 + # compute the average of the covariance matrices
1/M * lavInspect(path_fit,
"vcov.std.all")[c("Integrity_GB~TSM_GB",
"Integrity_CC~TSM_CC",
"Integrity_CB~TSM_CB"),
c("Integrity_GB~TSM_GB",
"Integrity_CC~TSM_CC",
"Integrity_CB~TSM_CB")]
}
# Compute the average of the estimates ###############
estimates_hyp2 <- colMeans(best_hyp2)
names(estimates_hyp2) <- c("Int_on_TSM_GB", "Int_on_TSM_CC", "Int_on_TSM_CB")
round(estimates_hyp2, 2)
## Int_on_TSM_GB Int_on_TSM_CC Int_on_TSM_CB
## -0.44 -0.41 -0.32
\usetikzlibrary{arrows,positioning}
\usetikzlibrary{calc}
\begin{tikzpicture}[auto,>=latex,align=center,
latent/.style={circle,draw, thick,inner sep=0pt,minimum size=10mm},
manifest/.style={rectangle,draw, thick,inner sep=2pt,minimum size=15mm},
manifestcov/.style={rectangle,draw, thick,inner sep=1pt,minimum size=8mm},
mean/.style={regular polygon,regular polygon sides=3,draw,thick,inner sep=0pt,minimum size=10mm},
paths/.style={->, thick, >=latex, font = \footnotesize\sffamily, fill = white},
variance/.style={<->, thin, >=latex, bend left=270, looseness=2.5, font = \footnotesize},
covariance/.style={<->, thin, >=latex, bend left=290, looseness=.5, font = \footnotesize},
dotsdots/.style={circle, draw, fill = black, minimum size = 1mm, maximum size = 1mm}
]
%% AVs
\node [manifest] (AV1) at (0, -2) {$int_{gb}$};
\node [manifest] (AV2) [right = of AV1] {$int_{cc}$};
\node [manifest] (AV3) [right = of AV2] {$int_{cb}$};
%% UVs
\node [manifest] (UV1) at (0, 2) {$tsm_{gb}$};
\node [manifest] (UV2) [right =of UV1] {$tsm_{cc}$};
\node [manifest] (UV3) [right =of UV2] {$tsm_{cb}$};
%% Paths
\draw [paths] (UV1) to node [midway]{-.44} (AV1);
\draw [paths] (UV2) to node [midway]{-.41} (AV2);
\draw [paths] (UV3) to node [midway]{-.32} (AV3);
%\draw [covariance] (UV2) to node {} (UV1);
\end{tikzpicture}
covbetween_hyp2 <- cov(best_hyp2) # between covariance matrix
covariance_hyp2 <- covwithin_hyp2 + (1+1/M)*covbetween_hyp2 # total variance
# Determine the effective and real sample sizes ######
samp_hyp2 <- nrow(data_scales_wide) # real sample size
nucom_hyp2 <- samp_hyp2-length(estimates_hyp2)
# corresponds to Equation (X) in Hoijtink, Gu, Mulder, & Rosseel (2019) ...
lam_hyp2 <- (1+1/M)*(1/length(estimates_hyp2))*
sum(diag(covbetween_hyp2 %*%
MASS::ginv(covariance_hyp2))) # ... (43)
nuold_hyp2 <- (M-1)/(lam_hyp2^2) # ... (44)
nuobs_hyp2 <- (nucom_hyp2+1)/(nucom_hyp2+3)*nucom_hyp2*(1-lam_hyp2) # ... (46)
nu_hyp2 <- nuold_hyp2*nuobs_hyp2/(nuold_hyp2+nuobs_hyp2) # ... (47)
fracmis_hyp2 <- (nu_hyp2+1)/(nu_hyp2+3)*lam_hyp2 + 2/(nu_hyp2+3) # ... (48)
neff_hyp2 <- samp_hyp2-samp_hyp2*fracmis_hyp2 # = 114 approx. 2/3* 270
# coerce `covariance` to a list
covariance_hyp2 <- list(covariance_hyp2)
set.seed(6346)
results_hyp2 <- bain(estimates_hyp2,
"Int_on_TSM_GB < 0 & Int_on_TSM_CC < 0 & Int_on_TSM_CB < 0;
Int_on_TSM_GB = 0 & Int_on_TSM_CC = 0 & Int_on_TSM_CB = 0",
Sigma = covariance_hyp2,
n = neff_hyp2,
group_parameters = 3,
joint_parameters = 0,
standardize = F)
print(results_hyp2)
## Bayesian informative hypothesis testing for an object of class numeric:
##
## Fit Com BF.u BF.c PMPa PMPb
## H1 1.000 0.161 6.225 10465021.465 1.000 0.862
## H2 0.000 1.281 0.000 0.000 0.000 0.000
## Hu 0.138
##
## Hypotheses:
## H1: Int_on_TSM_GB<0&Int_on_TSM_CC<0&Int_on_TSM_CB<0
## H2: Int_on_TSM_GB=0&Int_on_TSM_CC=0&Int_on_TSM_CB=0
##
## Note: BF.u denotes the Bayes factor of the hypothesis at hand versus the unconstrained hypothesis Hu. BF.c denotes the Bayes factor of the hypothesis at hand versus its complement.
results_hyp2$BFmatrix
## H1 H2
## H1 1.000000e+00 1.273079e+21
## H2 7.854974e-22 1.000000e+00
set.seed(6346)
results_hyp3 <- bain(estimates_hyp2,
"Int_on_TSM_GB < Int_on_TSM_CC < Int_on_TSM_CB;
(Int_on_TSM_GB, Int_on_TSM_CC) < Int_on_TSM_CB;
Int_on_TSM_GB = Int_on_TSM_CC = Int_on_TSM_CB",
Sigma = covariance_hyp2,
n = neff_hyp2,
group_parameters = 3,
joint_parameters = 0,
standardize = T)
print(results_hyp3)
## Bayesian informative hypothesis testing for an object of class numeric:
##
## Fit Com BF.u BF.c PMPa PMPb
## H1 0.547 0.174 3.142 5.733 0.127 0.122
## H2 0.815 0.343 2.375 8.431 0.096 0.092
## H3 10.003 0.518 19.310 19.310 0.778 0.748
## Hu 0.039
##
## Hypotheses:
## H1: Int_on_TSM_GB<Int_on_TSM_CC<Int_on_TSM_CB
## H2: (Int_on_TSM_GB,Int_on_TSM_CC)<Int_on_TSM_CB
## H3: Int_on_TSM_GB=Int_on_TSM_CC=Int_on_TSM_CB
##
## Note: BF.u denotes the Bayes factor of the hypothesis at hand versus the unconstrained hypothesis Hu. BF.c denotes the Bayes factor of the hypothesis at hand versus its complement.
results_hyp3$BFmatrix
## H1 H2 H3
## H1 1.0000000 1.322887 0.1627242
## H2 0.7559224 1.000000 0.1230069
## H3 6.1453658 8.129625 1.0000000
# Set up the matrices for the estimates ##############
mulest_hyp4 <- matrix(0,nrow=M,ncol=3) # setup of matrices to store multiple estimates
covwithin_hyp4 <- matrix(0,nrow=3,ncol=3) # and covariance matrices
# Estimate the coefficients for each data frame ######
for(i in 1:M) {
within_hyp4 <- lm(cbind( TSM_GB, TSM_CC, TSM_CB) ~ 1,
data=mice::complete(out,i)) # estimate the means of the three variables
mulest_hyp4[i,]<-coef(within_hyp4)[1:3] # store these means in the matrix `mulres`
covwithin_hyp4<-covwithin_hyp4 + 1/M * vcov(within_hyp4)[1:3,1:3] # compute the averages
}
# Compute the average of the estimates ###############
estimates_hyp4 <- colMeans(mulest_hyp4)
names(estimates_hyp4) <- c("TSM_GB","TSM_CC","TSM_CB")
covbetween_hyp4 <- cov(mulest_hyp4) # between covariance matrix
covariance_hyp4 <- covwithin_hyp4 + (1+1/M)*covbetween_hyp4 # total variance
# Determine the effective and real sample sizes ######
samp_hyp4 <- nrow(data_scales_wide) # real sample size
nucom_hyp4<-samp_hyp4-length(estimates_hyp4)
# corresponds to Equation (X) in Hoijtink, Gu, Mulder, & Rosseel (2019) ...
lam_hyp4 <- (1+1/M)*(1/length(estimates_hyp4))*
sum(diag(covbetween_hyp4 %*%
MASS::ginv(covariance_hyp4))) # ... (43)
nuold_hyp4<-(M-1)/(lam_hyp4^2) # ... (44)
nuobs_hyp4<-(nucom_hyp4+1)/(nucom_hyp4+3)*nucom_hyp4*(1-lam_hyp4) # ... (46)
nu_hyp4<- nuold_hyp4*nuobs_hyp4/(nuold_hyp4+nuobs_hyp4) # ... (47)
fracmis_hyp4 <- (nu_hyp4+1)/(nu_hyp4+3)*lam_hyp4 + 2/(nu_hyp4+3) # ... (48)
neff_hyp4<-samp_hyp4-samp_hyp4*fracmis_hyp4 # = 172 approx. 2/3* 270
# coerce `covariance` to a list
covariance_hyp4<-list(covariance_hyp4)
# Test the hypotheses with bain ######################
set.seed(6346)
results_hyp4 <- bain(estimates_hyp4,
"TSM_GB>TSM_CC>TSM_CB;
TSM_GB=TSM_CC=TSM_CB;
(TSM_GB,TSM_CC)>TSM_CB",
n = neff_hyp4, Sigma=covariance_hyp4,
group_parameters=3, joint_parameters = 0)
print(results_hyp4)
## Bayesian informative hypothesis testing for an object of class numeric:
##
## Fit Com BF.u BF.c PMPa PMPb
## H1 0.650 0.180 3.602 8.433 0.589 0.507
## H2 0.297 0.510 0.583 0.583 0.095 0.082
## H3 0.658 0.342 1.926 3.708 0.315 0.271
## Hu 0.141
##
## Hypotheses:
## H1: TSM_GB>TSM_CC>TSM_CB
## H2: TSM_GB=TSM_CC=TSM_CB
## H3: (TSM_GB,TSM_CC)>TSM_CB
##
## Note: BF.u denotes the Bayes factor of the hypothesis at hand versus the unconstrained hypothesis Hu. BF.c denotes the Bayes factor of the hypothesis at hand versus its complement.
results_hyp4$BFmatrix
## H1 H2 H3
## H1 1.0000000 6.181588 1.87024
## H2 0.1617707 1.000000 0.30255
## H3 0.5346909 3.305238 1.00000